﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Drawing;
using System.Drawing.Imaging;

namespace ImageProcessing
{
	public static class OrderStatFilters
	{


		internal static Bitmap Median(double[][][] adYIQ, int height, int width, int nWinLen, PixelFormat pixelformat)
		{
			int nWindowOffset = (nWinLen - 1) / 2;
			int nWinTotal = nWinLen * nWinLen;
			int k;
			int nRowWidth;
			int min;
			int half = nWinTotal / 2 + 1;
			int nWinI, nWinJ;
			double[][] adWindow = new double[nWinTotal][];
			double[] adPixel = new double[3];
			double[] adTemp;
			Bitmap bmOut = new Bitmap(width, height, pixelformat);

			#region Image pointer setup
			BitmapData srcBitData;
			Rectangle rectSrc = new Rectangle(0, 0, bmOut.Width, bmOut.Height);
			srcBitData = bmOut.LockBits(rectSrc, ImageLockMode.ReadOnly, bmOut.PixelFormat);
			int srcStride = srcBitData.Stride;

			IntPtr srcPtr = srcBitData.Scan0;
			int imageBPP = Image.GetPixelFormatSize(bmOut.PixelFormat);
			int srcBytePP = imageBPP / 8;
			#endregion Image pointer setup

			unsafe
			{
				byte* sptr = (byte*)(void*)srcPtr;

				//calculate values
				for (int j = 0; j < bmOut.Height; j++)
				{
					k = j * srcStride;
					for (int i = 0; i < bmOut.Width; i++)
					{
						if (j < nWindowOffset || j >= bmOut.Height - nWindowOffset || i < nWindowOffset || i >= bmOut.Width - nWindowOffset)
						{
							adPixel = ColorConversions.YIQConversion(adYIQ[i][j], false);
						}//end if

						else
						{
							nWinJ = j - nWindowOffset;

							//initialize whole window
							for (int n = 0; n < nWinLen; n++)
							{
								nRowWidth = n * nWinLen;
								nWinI = i - nWindowOffset;
								for (int m = 0; m < nWinLen; m++)
								{
									adWindow[nRowWidth + m] = adYIQ[nWinI][nWinJ];
									nWinI++;
								}
								nWinJ++;
							}


							//sort
							for (int m = 0; m < half; m++)
							{
								min = m;
								for (int n = m + 1; n < nWinTotal; n++)
								{
									if (adWindow[n][0] < adWindow[min][0])
									{
										min = n;
									}
								}
								//swap

								adTemp = adWindow[m];
								adWindow[m] = adWindow[min];
								adWindow[min] = adTemp;
								//only sort the first half because that's all we need for median

							}
							adPixel = ColorConversions.YIQConversion(adWindow[nWinTotal / 2], false);
						} //end else

						sptr[i * srcBytePP + k] = (byte)(adPixel[0] * 255.0);
						if (srcBytePP >= 3)
						{
							sptr[i * srcBytePP + k + 1] = (byte)(adPixel[1] * 255.0);
							sptr[i * srcBytePP + k + 2] = (byte)(adPixel[2] * 255.0);
							if (srcBytePP == 4)
								sptr[i * srcBytePP + k + 3] = 255;
						} //end if


					} //end for
				} //end for
			}
			bmOut.UnlockBits(srcBitData);
			return bmOut;

		}//end Median

		internal static Bitmap ArithmeticMean(double[][][] adYIQ, int height, int width, int nWinLen, PixelFormat pixelformat)
		{

			int nWindowOffset = (nWinLen - 1) / 2;
			int nWinTotal = nWinLen * nWinLen;
			int k;
			int half = nWinTotal / 2 + 1;
			int nWinI, nWinJ;
			double[][] adWindow = new double[nWinTotal][];
			double[] adPixel = new double[3];
			Bitmap bmOut = new Bitmap(width, height, pixelformat);

			#region Image pointer setup
			BitmapData srcBitData;
			Rectangle rectSrc = new Rectangle(0, 0, bmOut.Width, bmOut.Height);
			srcBitData = bmOut.LockBits(rectSrc, ImageLockMode.ReadOnly, bmOut.PixelFormat);
			int srcStride = srcBitData.Stride;

			IntPtr srcPtr = srcBitData.Scan0;
			int imageBPP = Image.GetPixelFormatSize(bmOut.PixelFormat);
			int srcBytePP = imageBPP / 8;
			#endregion Image pointer setup

			unsafe
			{
				byte* sptr = (byte*)(void*)srcPtr;

				//calculate values
				for (int j = 0; j < bmOut.Height; j++)
				{
					k = j * srcStride;
					for (int i = 0; i < bmOut.Width; i++)
					{
						if (j < nWindowOffset || j >= bmOut.Height - nWindowOffset || i < nWindowOffset || i >= bmOut.Width - nWindowOffset)
						{
							adPixel = ColorConversions.YIQConversion(adYIQ[i][j], false);
						}//end if

						else
						{
							nWinJ = j - nWindowOffset;

							//walk through window and sum intesity (Y) values
							adPixel[0] = 0;
							for (int n = 0; n < nWinLen; n++)
							{
								nWinI = i - nWindowOffset;
								for (int m = 0; m < nWinLen; m++)
								{
									adPixel[0] += adYIQ[nWinI][nWinJ][0];
									nWinI++;
								}
								nWinJ++;
							}
							adPixel[0] = adPixel[0] / nWinTotal; //calculate average 
							adPixel[1] = adYIQ[i][j][1]; //maintain same IQ values
							adPixel[2] = adYIQ[i][j][2];


							adPixel = ColorConversions.YIQConversion(adPixel, false);

						} //end else

						sptr[i * srcBytePP + k] = (byte)(adPixel[0] * 255.0);
						if (srcBytePP >= 3)
						{
							sptr[i * srcBytePP + k + 1] = (byte)(adPixel[1] * 255.0);
							sptr[i * srcBytePP + k + 2] = (byte)(adPixel[2] * 255.0);
							if (srcBytePP == 4)
								sptr[i * srcBytePP + k + 3] = 255;
						} //end if

					} //end for
				} //end for
			}
			bmOut.UnlockBits(srcBitData);
			return bmOut;
		}//end ArithmeticMean


		internal static Bitmap Max(double[][][] adYIQ, int height, int width, int nWinLen, PixelFormat pixelformat)
		{
			int nWindowOffset = (nWinLen - 1) / 2;
			int nWinTotal = nWinLen * nWinLen;
			int k;
			double[] adPixel = new double[3];
			int nWinI, nWinJ;
			Bitmap bmOut = new Bitmap(width, height, pixelformat);
			double max = double.MinValue;
			int[] maxpos = new int[2];

			#region Image pointer setup
			BitmapData srcBitData;
			Rectangle rectSrc = new Rectangle(0, 0, bmOut.Width, bmOut.Height);
			srcBitData = bmOut.LockBits(rectSrc, ImageLockMode.ReadOnly, bmOut.PixelFormat);
			int srcStride = srcBitData.Stride;

			IntPtr srcPtr = srcBitData.Scan0;
			int imageBPP = Image.GetPixelFormatSize(bmOut.PixelFormat);
			int srcBytePP = imageBPP / 8;
			#endregion Image pointer setup

			unsafe
			{
				byte* sptr = (byte*)(void*)srcPtr;

				//calculate values
				for (int j = 0; j < bmOut.Height; j++)
				{
					k = j * srcStride;
					for (int i = 0; i < bmOut.Width; i++)
					{
						if (j < nWindowOffset || j >= bmOut.Height - nWindowOffset || i < nWindowOffset || i >= bmOut.Width - nWindowOffset)
						{
							adPixel = ColorConversions.YIQConversion(adYIQ[i][j], false);
						}//end if

						else
						{
							nWinJ = j - nWindowOffset;
							//walk through window and sum intesity (Y) values
							adPixel[0] = 0;
							max = double.MinValue;
							for (int n = 0; n < nWinLen; n++)
							{
								nWinI = i - nWindowOffset;
								for (int m = 0; m < nWinLen; m++)
								{
									if (adYIQ[nWinI][nWinJ][0] > max)
									{
										max = adYIQ[nWinI][nWinJ][0];
										//maxpos[0] = nWinI;
										//maxpos[1] = nWinJ;
									} //end if 
								} //end for
								nWinJ++;
							} //end for

							adPixel[0] = max;
							adPixel[1] = adYIQ[i][j][1];
							adPixel[2] = adYIQ[i][j][2];

							adPixel = ColorConversions.YIQConversion(adPixel, false);
						} //end else

						sptr[i * srcBytePP + k] = (byte)(adPixel[0] * 255.0);
						if (srcBytePP >= 3)
						{
							sptr[i * srcBytePP + k + 1] = (byte)(adPixel[1] * 255.0);
							sptr[i * srcBytePP + k + 2] = (byte)(adPixel[2] * 255.0);
							if (srcBytePP == 4)
								sptr[i * srcBytePP + k + 3] = 255;
						} //end if

					} //end for
				} //end for
			}
			bmOut.UnlockBits(srcBitData);
			return bmOut;

		}//end Max

		internal static Bitmap GeometricMean(double[][][] adYIQ, int height, int width, int nWinLen, PixelFormat pixelformat)
		{

			int nWindowOffset = (nWinLen - 1) / 2;
			int nWinTotal = nWinLen * nWinLen;
			int k;
			double[] adPixel = new double[3];

			int nWinI, nWinJ;
			Bitmap bmOut = new Bitmap(width, height, pixelformat);
			int[] maxpos = new int[2];
			#region Image pointer setup
			BitmapData srcBitData;
			Rectangle rectSrc = new Rectangle(0, 0, bmOut.Width, bmOut.Height);
			srcBitData = bmOut.LockBits(rectSrc, ImageLockMode.ReadOnly, bmOut.PixelFormat);
			int srcStride = srcBitData.Stride;

			IntPtr srcPtr = srcBitData.Scan0;
			int imageBPP = Image.GetPixelFormatSize(bmOut.PixelFormat);
			int srcBytePP = imageBPP / 8;
			#endregion Image pointer setup

			unsafe
			{
				byte* sptr = (byte*)(void*)srcPtr;

				//calculate values
				for (int j = 0; j < bmOut.Height; j++)
				{
					k = j * srcStride;
					for (int i = 0; i < bmOut.Width; i++)
					{
						if (j < nWindowOffset || j >= bmOut.Height - nWindowOffset || i < nWindowOffset || i >= bmOut.Width - nWindowOffset)
						{
							adPixel = ColorConversions.YIQConversion(adYIQ[i][j], false);
						}//end if

						else
						{
							nWinJ = j - nWindowOffset;

							adPixel[0] = 1.0;

							//initialize whole window
							for (int n = 0; n < nWinLen; n++)
							{
								nWinI = i - nWindowOffset;
								for (int m = 0; m < nWinLen; m++)
								{
									adPixel[0] *= adYIQ[nWinI][nWinJ][0];

									nWinI++;
								}
								nWinJ++;
							}

							adPixel[0] = Math.Pow(adPixel[0], 1.0 / nWinTotal);
							adPixel[1] = adYIQ[i][j][1];
							adPixel[2] = adYIQ[i][j][2];

							adPixel = ColorConversions.YIQConversion(adPixel, false);
						} //end else

						sptr[i * srcBytePP + k] = (byte)(adPixel[0] * 255.0);
						if (srcBytePP >= 3)
						{
							sptr[i * srcBytePP + k + 1] = (byte)(adPixel[1] * 255.0);
							sptr[i * srcBytePP + k + 2] = (byte)(adPixel[2] * 255.0);
							if (srcBytePP == 4)
								sptr[i * srcBytePP + k + 3] = 255;
						} //end if

					} //end for
				} //end for
			}
			bmOut.UnlockBits(srcBitData);
			return bmOut;
		}//end GeometricMean

		internal static Bitmap HarmonicMean(double[][][] adYIQ, int height, int width, int nWinLen, PixelFormat pixelformat)
		{

			int nWindowOffset = (nWinLen - 1) / 2;
			int nWinTotal = nWinLen * nWinLen;
			int k;
			double[] adPixel = new double[3];

			int nWinI, nWinJ;
			Bitmap bmOut = new Bitmap(width, height, pixelformat);
			int[] maxpos = new int[2];
			#region Image pointer setup
			BitmapData srcBitData;
			Rectangle rectSrc = new Rectangle(0, 0, bmOut.Width, bmOut.Height);
			srcBitData = bmOut.LockBits(rectSrc, ImageLockMode.ReadOnly, bmOut.PixelFormat);
			int srcStride = srcBitData.Stride;

			IntPtr srcPtr = srcBitData.Scan0;
			int imageBPP = Image.GetPixelFormatSize(bmOut.PixelFormat);
			int srcBytePP = imageBPP / 8;
			#endregion Image pointer setup

			unsafe
			{
				byte* sptr = (byte*)(void*)srcPtr;

				//calculate values
				for (int j = 0; j < bmOut.Height; j++)
				{
					k = j * srcStride;
					for (int i = 0; i < bmOut.Width; i++)
					{
						if (j < nWindowOffset || j >= bmOut.Height - nWindowOffset || i < nWindowOffset || i >= bmOut.Width - nWindowOffset)
						{
							adPixel = ColorConversions.YIQConversion(adYIQ[i][j], false);
						}//end if

						else
						{
							nWinJ = j - nWindowOffset;

							adPixel[0] = 0.0;

							//initialize whole window
							for (int n = 0; n < nWinLen; n++)
							{
								nWinI = i - nWindowOffset;
								for (int m = 0; m < nWinLen; m++)
								{
									if (adYIQ[nWinI][nWinJ][0] > 0.0)
										adPixel[0] += 1.0/adYIQ[nWinI][nWinJ][0];

									nWinI++;
								}
								nWinJ++;
							}

							adPixel[0] = nWinTotal / adPixel[0];
							adPixel[1] = adYIQ[i][j][1];
							adPixel[2] = adYIQ[i][j][2];

							adPixel = ColorConversions.YIQConversion(adPixel, false);
						} //end else

						sptr[i * srcBytePP + k] = (byte)(adPixel[0] * 255.0);
						if (srcBytePP >= 3)
						{
							sptr[i * srcBytePP + k + 1] = (byte)(adPixel[1] * 255.0);
							sptr[i * srcBytePP + k + 2] = (byte)(adPixel[2] * 255.0);
							if (srcBytePP == 4)
								sptr[i * srcBytePP + k + 3] = 255;
						} //end if

					} //end for
				} //end for
			}
			bmOut.UnlockBits(srcBitData);
			return bmOut;
		}//end HarmonicMean

		internal static Bitmap Min(double[][][] adYIQ, int height, int width, int nWinLen, PixelFormat pixelformat)
		{
			int nWindowOffset = (nWinLen - 1) / 2;
			int nWinTotal = nWinLen * nWinLen;
			int k;
			double[] adPixel = new double[3];
			int nWinI, nWinJ;
			Bitmap bmOut = new Bitmap(width, height, pixelformat);
			double min = double.MaxValue;
			int[] minpos = new int[2];
			#region Image pointer setup
			BitmapData srcBitData;
			Rectangle rectSrc = new Rectangle(0, 0, bmOut.Width, bmOut.Height);
			srcBitData = bmOut.LockBits(rectSrc, ImageLockMode.ReadOnly, bmOut.PixelFormat);
			int srcStride = srcBitData.Stride;

			IntPtr srcPtr = srcBitData.Scan0;
			int imageBPP = Image.GetPixelFormatSize(bmOut.PixelFormat);
			int srcBytePP = imageBPP / 8;
			#endregion Image pointer setup

			unsafe
			{
				byte* sptr = (byte*)(void*)srcPtr;

				//calculate values
				for (int j = 0; j < bmOut.Height; j++)
				{
					k = j * srcStride;
					for (int i = 0; i < bmOut.Width; i++)
					{
						if (j < nWindowOffset || j >= bmOut.Height - nWindowOffset || i < nWindowOffset || i >= bmOut.Width - nWindowOffset)
						{
							adPixel = ColorConversions.YIQConversion(adYIQ[i][j], false);
						}//end if

						else
						{
							nWinJ = j - nWindowOffset;
							adPixel[0] = 0;
							min = double.MaxValue;
							for (int n = 0; n < nWinLen; n++)
							{
								nWinI = i - nWindowOffset;
								for (int m = 0; m < nWinLen; m++)
								{
									if (adYIQ[nWinI][nWinJ][0] < min)
									{
										min = adYIQ[nWinI][nWinJ][0];
										minpos[0] = nWinI;
										minpos[1] = nWinJ;
									} //end if 
								} //end for
								nWinJ++;
							} //end for

							adPixel[0] = adYIQ[minpos[0]][minpos[1]][0];
							adPixel[1] = adYIQ[minpos[0]][minpos[1]][1];
							adPixel[2] = adYIQ[minpos[0]][minpos[1]][2];

							adPixel = ColorConversions.YIQConversion(adPixel, false);
						} //end else

						sptr[i * srcBytePP + k] = (byte)(adPixel[0] * 255.0);
						if (srcBytePP >= 3)
						{
							sptr[i * srcBytePP + k + 1] = (byte)(adPixel[1] * 255.0);
							sptr[i * srcBytePP + k + 2] = (byte)(adPixel[2] * 255.0);
							if (srcBytePP == 4)
								sptr[i * srcBytePP + k + 3] = 255;
						} //end if

					} //end for
				} //end for
			}
			bmOut.UnlockBits(srcBitData);
			return bmOut;

		} //end Min


		internal static Bitmap MidPoint(double[][][] adYIQ, int height, int width, int nWinLen, PixelFormat pixelformat)
		{
			int nWindowOffset = (nWinLen - 1) / 2;
			int nWinTotal = nWinLen * nWinLen;
			int k;
			int nRowWidth;
			int nMinInd, nMaxInd;
			int nWinI, nWinJ;
			double[][] adWindow = new double[nWinTotal][];
			double[] adPixel = new double[3];
			Bitmap bmOut = new Bitmap(width, height, pixelformat);

			#region Image pointer setup
			BitmapData srcBitData;
			Rectangle rectSrc = new Rectangle(0, 0, bmOut.Width, bmOut.Height);
			srcBitData = bmOut.LockBits(rectSrc, ImageLockMode.ReadOnly, bmOut.PixelFormat);
			int srcStride = srcBitData.Stride;

			IntPtr srcPtr = srcBitData.Scan0;
			int imageBPP = Image.GetPixelFormatSize(bmOut.PixelFormat);
			int srcBytePP = imageBPP / 8;
			#endregion Image pointer setup

			unsafe
			{
				byte* sptr = (byte*)(void*)srcPtr;

				//calculate values
				for (int j = 0; j < bmOut.Height; j++)
				{
					k = j * srcStride;
					for (int i = 0; i < bmOut.Width; i++)
					{
						if (j < nWindowOffset || j >= bmOut.Height - nWindowOffset || i < nWindowOffset || i >= bmOut.Width - nWindowOffset)
						{
							adPixel = ColorConversions.YIQConversion(adYIQ[i][j], false);
						}//end if

						else
						{
							nWinJ = j - nWindowOffset;

							//initialize whole window
							for (int n = 0; n < nWinLen; n++)
							{
								nRowWidth = n * nWinLen;
								nWinI = i - nWindowOffset;
								for (int m = 0; m < nWinLen; m++)
								{
									adWindow[nRowWidth + m] = adYIQ[nWinI][nWinJ];
									nWinI++;
								}
								nWinJ++;
							}

							//find min and max
							nMinInd = nMaxInd = 0;
							for (int m = 1; m < nWinTotal; m++)
							{
								if (adWindow[m][0] < adWindow[nMinInd][0])
									nMinInd = m;
								if (adWindow[m][0] > adWindow[nMaxInd][0])
									nMaxInd = m;
							} //end for

							
							adPixel[0] = (adWindow[nMinInd][0] + adWindow[nMaxInd][0]) / 2.0;
							adPixel[1] = adYIQ[i][j][1];
							adPixel[2] = adYIQ[i][j][2];
							adPixel = ColorConversions.YIQConversion(adPixel, false);
						} //end else

						sptr[i * srcBytePP + k] = (byte)(adPixel[0] * 255.0);
						if (srcBytePP >= 3)
						{
							sptr[i * srcBytePP + k + 1] = (byte)(adPixel[1] * 255.0);
							sptr[i * srcBytePP + k + 2] = (byte)(adPixel[2] * 255.0);
							if (srcBytePP == 4)
								sptr[i * srcBytePP + k + 3] = 255;
						} //end if

					} //end for
				} //end for
			}
			bmOut.UnlockBits(srcBitData);
			return bmOut;

		} //end Midpoint

		internal static Bitmap AlphaTrimMean(double[][][] adYIQ, int height, int width, int nWinLen, PixelFormat pixelformat, int alpha)
		{
			int nWindowOffset = (nWinLen - 1) / 2;
			int nWinTotal = nWinLen * nWinLen;
			int k;
			int nRowWidth;
			int nMinInd, nUpperInd;
			int nHalf;
			int nWinI, nWinJ;
			double[][] adWindow = new double[nWinTotal][];
			double[] adPixel = new double[3];
			double[] adTemp;

			if (alpha > nWinTotal)
			{
				alpha = nWinTotal - 1;
			} //end if
			Bitmap bmOut = new Bitmap(width, height, pixelformat);

			#region Image pointer setup
			BitmapData srcBitData;
			Rectangle rectSrc = new Rectangle(0, 0, bmOut.Width, bmOut.Height);
			srcBitData = bmOut.LockBits(rectSrc, ImageLockMode.ReadOnly, bmOut.PixelFormat);
			int srcStride = srcBitData.Stride;

			IntPtr srcPtr = srcBitData.Scan0;
			int imageBPP = Image.GetPixelFormatSize(bmOut.PixelFormat);
			int srcBytePP = imageBPP / 8;
			#endregion Image pointer setup

			unsafe
			{
				byte* sptr = (byte*)(void*)srcPtr;

				//calculate values
				for (int j = 0; j < bmOut.Height; j++)
				{
					k = j * srcStride;
					for (int i = 0; i < bmOut.Width; i++)
					{
						if (j < nWindowOffset || j >= bmOut.Height - nWindowOffset || i < nWindowOffset || i >= bmOut.Width - nWindowOffset)
						{
							adPixel = ColorConversions.YIQConversion(adYIQ[i][j], false);
						}//end if

						else
						{
							nWinJ = j - nWindowOffset;
							
							//initialize window
							for (int n = 0; n < nWinLen; n++)
							{
								nRowWidth = n * nWinLen;
								nWinI = i - nWindowOffset;
								for (int m = 0; m < nWinLen; m++)
								{
									adWindow[nRowWidth + m] = adYIQ[nWinI][nWinJ];
									nWinI++;
								}
								nWinJ++;
							}

							//sort
							for (int m = 0; m < nWinTotal; m++)
							{
								nMinInd = m;
								for (int n = m + 1; n < nWinTotal; n++)
								{
									if (adWindow[n][0] < adWindow[nMinInd][0])
									{
										nMinInd = n;
									}
								}

								//swap
								if (nMinInd != m)
								{
									adTemp = adWindow[m];
									adWindow[m] = adWindow[nMinInd];
									adWindow[nMinInd] = adTemp;
								} //end if 
							}

							adPixel[0] = 0;
							nHalf = nWinTotal / 2;
							nUpperInd = nHalf + 1;
							for (int m = 0; m < nWinTotal - alpha; m++)
							{
								if (m % 2 == 0)
								{
									adPixel[0] += adWindow[nHalf][0];
									nHalf--;
								} //end if
								else
								{
									adPixel[0] += adWindow[nUpperInd][0];
									nUpperInd++;
								} //end if
							} //end for

							adPixel[0] = adPixel[0] / (nWinTotal - alpha);
							adPixel[1] = adYIQ[i][j][1];
							adPixel[2] = adYIQ[i][j][2];
							adPixel = ColorConversions.YIQConversion(adPixel, false);
						} //end else

						sptr[i * srcBytePP + k] = (byte)(adPixel[0] * 255.0);
						if (srcBytePP >= 3)
						{
							sptr[i * srcBytePP + k + 1] = (byte)(adPixel[1] * 255.0);
							sptr[i * srcBytePP + k + 2] = (byte)(adPixel[2] * 255.0);
							if (srcBytePP == 4)
								sptr[i * srcBytePP + k + 3] = 255;
						} //end if

					} //end for
				} //end for
			}
			bmOut.UnlockBits(srcBitData);
			return bmOut;

		} //end AlphaTrimMean

		internal static Bitmap Contraharmonic(double[][][] adYIQ, int height, int width, int nWinLen, PixelFormat pixelformat, int q)
		{

			int nWindowOffset = (nWinLen - 1) / 2;
			int nWinTotal = nWinLen * nWinLen;
			int k;
			int half = nWinTotal / 2 + 1;
			int nWinI, nWinJ;
			double[][] adWindow = new double[nWinTotal][];
			double[] adPixel = new double[3];
			Bitmap bmOut = new Bitmap(width, height, pixelformat);

			#region Image pointer setup
			BitmapData srcBitData;
			Rectangle rectSrc = new Rectangle(0, 0, bmOut.Width, bmOut.Height);
			srcBitData = bmOut.LockBits(rectSrc, ImageLockMode.ReadOnly, bmOut.PixelFormat);
			int srcStride = srcBitData.Stride;

			IntPtr srcPtr = srcBitData.Scan0;
			int imageBPP = Image.GetPixelFormatSize(bmOut.PixelFormat);
			int srcBytePP = imageBPP / 8;
			#endregion Image pointer setup

			unsafe
			{
				byte* sptr = (byte*)(void*)srcPtr;

				//calculate values
				for (int j = 0; j < bmOut.Height; j++)
				{
					k = j * srcStride;
					for (int i = 0; i < bmOut.Width; i++)
					{
						if (j < nWindowOffset || j >= bmOut.Height - nWindowOffset || i < nWindowOffset || i >= bmOut.Width - nWindowOffset)
						{
							adPixel = ColorConversions.YIQConversion(adYIQ[i][j], false);
						}//end if

						else
						{
							nWinJ = j - nWindowOffset;

							//walk through window and sum intesity (Y) values
							adPixel[0] = 0;
							adPixel[1] = 0;
							for (int n = 0; n < nWinLen; n++)
							{
								nWinI = i - nWindowOffset;
								for (int m = 0; m < nWinLen; m++)
								{
									adPixel[0] += Math.Pow(adYIQ[nWinI][nWinJ][0], q + 1);
									adPixel[1] += Math.Pow(adYIQ[nWinI][nWinJ][0],q);
									nWinI++;
								}
								nWinJ++;
							}
							adPixel[0] = adPixel[0] / adPixel[1];
							adPixel[1] = adYIQ[i][j][1]; //maintain same IQ values
							adPixel[2] = adYIQ[i][j][2];


							adPixel = ColorConversions.YIQConversion(adPixel, false);

						} //end else

						sptr[i * srcBytePP + k] = (byte)(adPixel[0] * 255.0);
						if (srcBytePP >= 3)
						{
							sptr[i * srcBytePP + k + 1] = (byte)(adPixel[1] * 255.0);
							sptr[i * srcBytePP + k + 2] = (byte)(adPixel[2] * 255.0);
							if (srcBytePP == 4)
								sptr[i * srcBytePP + k + 3] = 255;
						} //end if

					} //end for
				} //end for
			}
			bmOut.UnlockBits(srcBitData);
			return bmOut;
		}//end Contraharmonic

	}
}
